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Analyzing  and  Exploiting  Numerical  Characteristics  of 

Zone  Fire  Models 


Glenn  P.  Forney*  William  F.  Moss^ 


Abstract 

In  order  to  design  robust  and  stable  zone  fire  modeling  algorithms,  the  nu- 
merical properties  of  computer  arithmetic  and  the  modeling  differential  equa- 
tions must  be  understood.  This  report  examines  some  of  these  properties  and 
provides  tools  for  their  analysis.  Many  sets  of  differential  equations  for  zone  fire 
modeling  can  be  derived  using  the  laws  of  conservation  of  mass  and  energy.  A 
comparison  between  various  possible  formulations  is  made  in  terms  of  numer- 
ical properties.  One  property  that  many  formulations  possess  is  the  presence 
of  multiple  time  scales.  Pressures  equilibrate  much  faster  than  other  quanti- 
ties such  as  density  and  temperature.  Numerically,  this  property  is  known  as 
stiffness.  Stiffness,  in  the  context  of  fire  modeling,  and  numerical  methods  for 
handling  it  are  discussed. 


1 Introduction 

1.1  Background 

An  understanding  of  the  numerical  properties  of  computer  arithmetic  and  the  dif- 
ferential equations  used  in  modeling  is  required  for  the  design  of  robust  and  stable 
zone  fire  modeling  algorithms.  The  basic  premise  used  to  formulate  a zone  fire 
model  is  that  an  enclosure  can  be  divided  into  a number  of  regions  or  zones  each 
with  approximately  uniform  conditions.  These  zones  interact  by  exchanging  mass 
and  energy [1].  Mass  and  energy  conservation  along  with  expressions  relating  mass, 
density,  volume,  internal  energy,  temperature  and  pressure  can  be  used  to  show 
that  many  formulations  exist  for  tracking  conditions  in  zones.  These  formulations 
are  equivalent  in  the  sense  that  one  formulation  may  be  converted  to  another  using 
physical  laws  such  as  the  ideal  gas  law  or  definitions  such  as  that  for  density  or  inter- 
nal energy.  Computationally,  zone  fire  modeling  is  challenging  due  to  the  numerical 

’Building  and  Fire  Research  Laboratory,  National  Institute  of  Standards  and  Technology, 
Gaithersburg,  MD  20899,  U.S.A.  (gpfrnCcfr6.cfr.nist.gov). 

'Department  of  Mathematical  Sciences,  Clemson  University,  Clemson,  SC  29634-1907,  U.S.A. 
(bmossCmath . clemson . edu) . 


1 


characteristics  of  the  basic  conservation  equations  used  to  simulate  mass  and  energy 
exchange  between  various  zones.  The  purpose  of  this  report  then  is  to  provide  a 
numerical  foundation  for  the  design  of  fire  modeling  algorithms.  It  is  important  to 
understand  when  differences  in  these  algorithms  are  numerically  significant  so  that 
the  best  possible  fire  modeling  algorithms  can  be  designed  and  implemented. 

1.2  Overview 

Fire  modeling  algorithms  need  to  be  designed  with  the  limitations  of  floating  point 
arithmetic  in  mind.  The  types  of  error  encountered  in  a computational  algorithm  are 
described.  Basic  tools  from  numerical  analysis  are  presented  in  order  to  understand 
how  errors  introduced  in  the  modeling  process  propagate  or  grow. 

Several  zone  fire  modeling  formulations  are  derived  using  conservation  of  mass 
and  energy.  The  advantages  and  disadvantages  of  each  formulation  from  a numerical 
perspective  are  indicated.  Approximate  formulations  based  on  the  rapid  equilibra- 
tion of  pressure  are  discussed.  Other  approximate  formulations  based  on  a simplified 
treatment  of  the  lower  layer  are  also  discussed. 

The  numerical  characteristics  of  a few  physical  models  for  exchanging  mass  and 
energy  between  zones  where  interesting  numerical  problems  have  been  observed 
are  discussed.  The  behavior  of  the  differential  equations  treated  as  a system  is 
also  discussed.  The  fact  that  pressure  in  a room  equilibrates  rapidly  in  response  to 
changing  fire  size,  vent  area,  and  layer  height  causes  difficulties  in  the  solution  of  the 
differential  equations.  This  phenomenon  is  called  stiffness.  Algorithms  for  solving 
stiff  systems  of  ordinary  differential  equations  are  outlined. 

Floating  point  numbers,  model  numbers  and  machine  or  computer  numbers  are 
all  terms  that  refer  to  those  numbers  that  are  exactly  representable  on  the  com- 
puter. The  appendix  presents  a model  of  floating  point  arithmetic  which  details 
how  numbers  are  represented  in  computers  and  how  operations  with  these  numbers 
behave. 

2 Sources  and  Analysis  of  Numerical  Error 

The  introduction  of  error  into  a fire  modeling  algorithm  is  inevitable;  however,  the 
algorithm  should  be  designed  to  minimize  the  impact  of  these  errors.  Some  common 
sources  of  error  are  physical  models  which  do  not  completely  describe  the  phenomena 
of  interest,  imprecise  data,  limitations  of  the  algorithms  used  to  solve  the  modeling 
equations,  and  numerical  errors  introduced  when  these  algorithms  are  implemented 
on  a computer.  This  report  is  mainly  concerned  with  the  last  two  sources  of  error. 

2.1  Sources  of  Numerical  Error 

Numerical  error  can  be  divided  into  three  categories:  roundoff,  truncation  and  dis- 
cretization error.  Roundoff  error  occurs  because  computers  represent  real  numbers 
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using  a finite  number  of  digits.  The  best  that  can  be  expected  is  that  the  computed 
result  of  an  operation  is  the  nearest  floating  point  number  to  the  true  result.  The 
appendix  details  the  properties  of  computer  arithmetic  as  implemented  on  most 
computers  used  today.  These  properties  can  be  summarized  by  the  statement  that 
“subtractions  can  cause  significant  loss  of  accuracy  in  a numerical  computation  and 
this  loss  of  accuracy  can  be  greatly  amplified  in  subsequent  calculations,  a condi- 
tion termed  catastrophic  cancellation.”  Truncation  error  occurs  when  an  infinite 
process  is  replaced  by  a finite  one.  This  can  happen,  for  example,  when  an  infinite 
series  is  truncated  after  a finite  number  of  terms  or  when  an  iteration  is  terminated 
after  a convergence  criterion  has  been  satisfied.  Discretization  error  occurs  when 
a continuous  process  such  as  a derivative  is  approximated  by  a discrete  analog  such 
as  a divided  difference. 

Consider  the  problem  of  finding  the  solution  at  f > 0 to  the  scalar  initial  value 
problem 

fr  = 

y(0)  = h . 

Denote  this  solution  by  y{t).  A discrete  analog  of  problem  (1)  is  obtained  by  re- 
placing the  derivative  by  a divided  difference.  The  simplest  example  is  the  Euler 
method.  The  advantages  and  disadvantages  discussed  below  of  using  the  explicit  and 
implicit  Euler  methods  are  representative  of  more  complicated  methods  for  solving 
differential  equations.  These  alternate  methods  differ  in  the  number  of  terms  in  the 
Taylor  expansion  that  are  considered. 

Choose  a positive  integer  n (the  number  of  steps)  and  define  r,-  = ih,  z = 1, . . . , n, 
where  the  stepsize  h = t/n.  Then  = t and  y{i)  is  approximated  by  which  is 
defined  by  the  recurrence 

= fin,  yi),  z = 1, . . . , n - 1 

yo  = b . 

The  quantity  yiTi)  — y{  is  usually  called  the  global  error.  Disregarding  roundoff  error, 
the  final  global  error,  y{t)  — z/„,  can  be  shown  to  be  proportional  to  the  stepsize  h 
under  fairly  general  assumptions  about  / (see  [2]). 

A problem  such  as  (1)  is  called  stable  or  well-conditioned  if  small  changes 
in  the  data  or  input  parameters  produce  small  changes  in  the  solution;  otherwise, 
the  problem  is  called  unstable  or  ill-conditioned.  It  is  important  to  note  that 
this  concept  has  nothing  to  do  with  numerical  methods  for  solving  the  problem.  If 
fi'’’,y)  = o.y  ill  (1)  with  a > 0,  then  this  problem  is  unstable  or  ill-conditioned.  This 
can  be  seen  by  examining  the  family  of  solutions  to  the  differential  equation  that 
pass  close  to  the  initial  point  (0,  h).  As  illustrated  in  Figure  1,  two  of  these  solutions 
can  be  close  at  t = 0,  but  widely  separated  at  t > 0.  Similarly,  if  /(r,  y)  = ay  in  (1) 
with  a < 0,  then  this  problem  is  stable  or  weU-conditioned. 
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Figure  1:  A Graphical  Example  of  Instability 


The  inherent  error  of  a problem  is  the  maximum  error  that  could  occur  in  the 
solution  if  the  data  were  perturbed  by  a relative  amount  no  greater  than  machine 
epsilon  where  machine  epsilon  is  the  smallest  positive  number  added  to  one  that  gives 
a result  different  than  one.  This  and  other  floating  point  properties  are  discussed 
in  more  detail  in  the  appendix.  In  the  simplest  case  where  f{T,y)  = ay  in  (1),  the 
inherent  error  would  arise  from  perturbations  in  a and  b.  An  algorithm  for  solving 
a problem  is  called  numerically  stable  if  the  error  in  the  numerical  solution  is  the 
same  order  of  magnitude  as  the  inherent  error;  otherwise,  it  is  called  numerically 
unstable.  Roughly  speaking,  an  algorithm  is  numerically  unstable  if  small  errors 
made  at  one  stage  are  amplified  in  subsequent  stages  and  seriously  degrade  the 
accuracy  of  the  overall  computation. 

It  is  important  to  distinguish  between  iU-conditioned  problems  and  numerically 
unstable  algorithms.  An  iU- conditioned  problem  cannot  be  made  well-conditioned 
by  simply  using  stable  algorithms.  AU  algorithms  will  have  difficulties  for  such 
problems.  In  extreme  cases,  when  increasing  the  precision  of  the  computation  does 
not  help,  it  may  be  necessary  to  modify  the  problem. 

Consider  the  Euler  method  applied  to  (1)  with  /(t,  y)  = —ay  and  a > 0.  The 
analytic  solution  is  y(t)  = yoexp{—at)  which  decays  with  increasing  t.  The  discrete 
solution  using  the  explicit  Euler  method  and  assuming  exact  arithmetic  is  t/,-  = 
yi-i(l  - ah)  = yo{l  - ah)'.  The  growth  of  the  global  error,  e,-  :=  |y(r,)  - y,|,  will 
determine  if  the  Euler  method  is  numerically  stable.  The  global  error  at  step  i 
can  be  expressed  in  terms  of  the  global  error  at  the  previous  step  to  be 

e,+i  = |y(r,+i)  - j/i+il 

= - (1  - (^h))y{Ti)  + (1  - ah){y{Ti)  - y,)| 

< yo^-^  - ah\ei  . 

The  expression  |1  — a/i|  is  called  the  amplification  factor.  If  h > 2/a,  then  |1  — a/i|  > 1 
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and  both  e,  and  grow  exponentially  with  the  index  i.  In  this  case  the  Euler  method 
will  be  numerically  unstable  (even  on  a computer  that  uses  an  infinite  number  of 
digits  in  its  arithmetic). 

If  a > 0,  the  best  option  is  to  change  to  an  implicit  method  such  as  the  backward 
Euler  method  given  by 


Vi+i  - Vi 
h 

yo 


f{ri+i,yi+i),  i = I 

b . 


The  error  bound  for  the  backward  Euler  method  is 


C:+l 


< 


|1  + ah\ 


Here  the  amplification  factor  is  1/(1  + ah).  For  aU  positive  h,  the  amplification 
factor  is  less  than  or  equal  to  one  and  the  backward  Euler  method  is  numerically 
stable. 

This  example  says  that  an  unstable  algorithm  can  give  the  wrong  result  for 
a stable  or  well-conditioned  problem,  whereas  a stable  algorithm  will  give  correct 
results.  Further,  ill-conditioned  problems  cannot  be  made  weU-conditioned  by  using 
stable  algorithms.  Strategies  for  improving  the  solution  of  ill-conditioned  problems 
must  center  on  methods  for  changing  the  problem  to  make  it  more  stable. 


2.2  Analyzing  Numerical  Error 

Floating  point  error  bounds,  such  as  those  presented  in  the  appendix,  derived  using 
fundamental  properties  of  computer  arithmetic,  are  difficult  to  apply  to  complicated 
problems.  Simpler  tools,  however,  can  be  derived  for  analyzing  error  propagation 
properties  of  differentiable  functions.  Suppose  that  a fire  algorithm  can  be  expressed 
as  a function  f{x).  As  an  example,  the  input  x could  represent  pressure  and  the 
value  of  / could  represent  enthalpy  flow  rate  through  a vent.  If  x is  perturbed  by 
an  amount  h,  how  is  f{x)  affected?  This  question  can  be  addressed  using  condition 
numbers.  Suppose  / has  two  continuous  derivatives.  From  Taylor’s  Theorem,  it 
follows  that 

Af  :=  fix  + h)-  fix)  = f'ix)h  -b  0(h2)  (2) 

where  Oih?)  denotes  a term  that  is  bounded  by  a constant  times  h^.  The  absolute 
condition  number  of  /,  denoted  Ca(/),  is  defined  to  be  the  coefficient  of  h in  an  ex- 
pansion of  Af  in  powers  of  h.  To  a first  order  approximation,  the  absolute  condition 
number  relates  the  absolute  changes  or  error  in  x with  the  absolute  changes  in  / 
according  to 


A/  = Caif)h. 
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Consequently 


Caif)  = fix)  . 

Similarly,  the  relative  condition  number  of  f,  denoted  Cr(/),  is  defined  to  be  the 
coefficient  of  h/x  in  an  expansion  of  A///  in  powers  of  h/x.  From  (2)  it  follows 
that 


Crif) 


xfjx) 

fi^) 


Similarly,  the  relative  condition  number  relates  the  relative  error  of  x with  the 
relative  error  of  / according  to 


— = 


The  absolute  and  relative  condition  numbers  Ca  and  Cr  can  be  used  to  analyze 
numerical  properties  of  fire  algorithms.  These  condition  numbers  give  a quantitative 
measure  of  how  input  errors  are  magnified  by  a problem.  For  example,  consider  the 
expression 


9vent  — K\/Pi  — P2 


for  enthalpy  flow  rate  through  a vent  where  K = CpC vent ^ vent y/‘^PventTvent » Cp  is  the 
constant  pressure  specific  heat,  Avent  is  the  area  of  the  vent,  Cvent  is  the  vent  flow 
coefficient,  and  Pvent  a-nd  Tvent  are  the  density  and  temperature  of  the  gas  flowing 
through  the  vent.  The  pressure  difference  Pi  — P2  across  the  vent  drives  the  flow  of 
mass  through  the  vent.  The  relative  and  absolute  condition  numbers  for  ^vent  with 
respect  to  pressure  Pi  are 


^ai^vent)  — 

^r(9vent)  — 

Both  condition  numbers  show  that  problems  can  occur  when  Pi  and  P2  are  close 
and  in  the  case  of  the  relative  condition  number,  these  problems  are  independent  of 
K. 

Condition  numbers  can  be  generalized  to  predict  error  propagation  properties  of 
vector- valued  functions,  for  example  see  [3].  Table  1 gives  the  absolute  and  relative 
condition  numbers  for  a few  elementary  functions. 


9 vent 

Pi 

Pi  -P2' 
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Table  1:  Absolute  and  Relative  Condition  Numbers  for  a Few  Elementary  Functions 


/ 

Ca(/) 

Crif) 

fix) 

fix) 

f(^) 

X a 

1 

X 

x+a 

ax 

a 

1 

sin(a;) 

cos(x) 

X cot(x) 

a;“ 

a 

eax 

ax 

3 Zone  Fire  Modeling  Formulations 

The  zone  fire  models  presented  here  take  the  mathematical  form  of  an  initial  value 
problem  for  a system  of  differential  equations.  These  equations  are  derived  using  the 
conservation  of  mass  or  continuity  equation,  the  conservation  of  energy  or  the  first 
law  of  thermodynamics,  the  ideal  gas  law,  and  definitions  of  density  and  internal 
energy  (for  example,  see  [4]).  The  conservation  of  momentum  is  ignored.  These 
conservation  laws  are  invoked  for  each  zone  or  control  volume.  A zone  may  consist 
of  a number  of  interior  regions  (usually  an  upper  and  a lower  gas  layer),  and  a number 
of  wall  segments.  The  basic  assumption  of  a zone  fire  model  is  that  properties  such  as 
temperatures  can  be  uniformly  approximated  throughout  the  zone.  It  is  remarkable 
that  this  assumption  seems  to  hold  for  as  few  as  two  gas  layers. 

Many  differential  equation  formulations  based  upon  these  assumptions  can  be 
derived.  One  formulation  can  be  converted  into  another  using  definitions  of  density, 
internal  energy  and  the  ideal  gas  law.  Though  equivalent  analytically,  these  formu- 
lations differ  in  their  numerical  properties.  One  property  that  many  share  is  the 
presence  of  multiple  time  scales.  Physically,  the  pressure  in  a compartment  equili- 
brates much  quicker  than  densities  and  temperatures.  Numerically,  this  property  is 
known  as  stiffness  and  requires  the  use  of  special  differential  equation  solvers.  The 
physical  origins  and  numerical  consequences  of  stiffness  are  discussed  in  more  detail 
in  section  4.2. 

Each  differential  formulation  can  be  expressed  in  terms  of  mass  and  enthalpy  flow 
rates.  These  flow  rates  represent  the  exchange  of  mass  and/or  energy  between  zones 
due  to  physical  phenomena  or  sub-models  such  as  fire  plumes,  natural  and  forced 
vents,  convective  and  radiative  heat  transfer  etc.  For  example,  a vent  exchanges 
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mass  and  energy  between  zones  in  connected  rooms,  a fire  plume  typically  adds  heat 
to  the  upper  layer  and  transfers  entrained  mass  and  energy  from  the  lower  to  the 
upper  layer,  and  convection  transfers  energy  from  the  gas  layers  to  the  surrounding 
walls. 

Using  the  formalism  developed  by  Cooper  for  CCFM. VENTS  [4]  the  mass  flow 
rate  to  the  upper  and  lower  layers  is  denoted  mu  and  mi  and  the  enthalpy  flow 
rate  to  the  upper  and  lower  layers  is  denoted  qu  and  qi.  It  is  tacitly  assumed  that 
these  flow  rates  may  be  computed  in  terms  of  zone  properties  such  as  temperatures, 
densities,  etc.  These  rates  represent  the  net  sum  of  all  possible  sources  of  mass  and 
energy  due  to  phenomena  like  those  listed  above.  The  numerical  characteristics  of 
the  differential  equation  formulations  are  easier  to  identify  if  the  underlying  physical 
phenomena  are  decoupled  in  this  way. 

Many  approximations  are  obviously  necessary  when  developing  physical  sub- 
models for  the  mass  and  enthalpy  flow  rate  terms.  For  example,  most  fire  models 
assume  that  1)  the  specific  heat  terms  Cp  and  Cy  are  constant  even  though  they 
depend  upon  temperature,  2)  hydrostatic  terms  can  be  ignored  in  the  equation  of 
state  (the  ideal  gas  law)  relating  density  of  a layer  with  its  temperature.  We  wish  to 
distinguish  between  various  formulations  according  to  whether  they  are  mathemati- 
cally equivalent  to  the  conservation  laws  of  mass  and  energy.  A formulation  which  is 
equivalent  to  the  conservation  laws  will  be  denoted  conservative  otherwise  it  will 
be  identified  as  approximate.  Conservative  formulations  in  this  sense  are  not  nec- 
essarily better  than  approximate  ones.  The  next  two  sections  discuss  formulations 
which  are  conservative  and  approximate.  Again,  two  conservative  formulations  that 
are  equivalent  mathematically  need  not  be  equivalent  numerically. 

3.1  Conservative  Formulations 

A compartment  can  be  divided  into  two  control  volumes,  an  upper  layer  of  hot  gases 
and  smoke  and  a lower  layer  of  air  as  illustrated  in  Figure  2.  The  gas  in  each  layer 
has  attributes  of  mass,  internal  energy,  density,  temperature,  and  volume  denoted 
respectively  by  m,-,  Ei,  pi,  Ti,  and  Vi  where  i = L for  the  lower  layer  and  i = U 
for  the  upper  layer.  The  compartment  as  a whole  has  the  attribute  of  pressure  P. 
These  eleven  variables  are  related  by  means  of  the  following  seven^  constraints 


Pi  = 

Y (density) 

(3) 

Ei  = 

CymiTi  (internal  energy) 

(4) 

P = 

RpiTi  (ideal  gas  law) 

(5) 

V = 

Fl  + Vu  (total  volume)  . 

(6) 

^We  get  seven  by  counting  density, 

internal  energy  and  the  ideal 

gas  law  twice  (once  for  each 

layer). 
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Components  of  mass  and 
enthdpy  entering  or  leaving  a 
zone 


Figure  2:  Two  Layer  Zone  Model  Setup 


The  specific  heats  at  constant  volume  and  at  constant  pressure  Cy  and  Cp,  the  uni- 
versal gas  constant,  i?,  and  the  ratio  of  specific  heats,  7,  are  related  by 


7 = 
R = 


For  air,  Cp  « 1000  kJ/kg  K and  7 = 1.4. 

The  differential  equations  for  mass  in  each  layer  are  trivially 


dmi 

dt 

dmu 

dt 


mi 


mu  . 


The  first  law  of  thermodynamics  states  that  the  rate  of  increase  of  layer  internal 
energy  plus  the  rate  at  which  the  layer  does  work  by  expansion  is  equal  to  the  rate 
at  which  enthalpy  is  added  to  the  gas.  In  differential  equation  form  this  is 


internal  energy  work 


dt  dt 


enthalpy 


(7) 


A differential  equation  for  pressure  can  be  derived  by  adding  the  upper  and  lower 


9 


layer  versions  of  equation  (7),  noting  that 


dVj, 

dt 


dVr, 

dt 


and 


to  obtain 


dEi  di^c.ijTTiiTi)  Cfj  d f p-y  \ 

dt  dt  R dt  * 


dP  7 - 1..  , . 

lu  = • 


(8) 


(9) 


Differential  equations  for  the  layer  volumes  can  be  obtained  by  substituting  equation 
(8)  into  equation  (7)  to  obtain 

f = ^((7  - IW.  - Kf ) . (10) 

Equation  (7)  can  be  rewritten  to  eliminate  the  ^ term  to  obtain 

dEi  1 jrdP. 

Differential  equations  for  the  densities  can  be  derived  by  applying  the  quotient  rule 
to  ^ ^ (^)  using  equation  (10)  to  eliminate  ^ to  obtain 


dpi 

dt 


-1  ^ 


(11) 


Differential  equations  for  temperatures  can  be  derived  by  applying  the  quotient  rule 
to  ^ ^ using  equation  (11)  to  eliminate  ^ to  obtain 


dTi  1 //•  • T.,.rdP. 


(12) 


Differential  equations  for  each  of  the  eleven  variables  are  summarized  in  Table  2. 
Notice  that  a ^ term  occurs  in  all  but  the  mass  equations.  For  many  fire  scenarios 
this  term  can  be  set  to  zero.  Section  3.3  discusses  approximations  to  the  zone  fire 
modeling  equations  derived  by  dropping  the  pressure  transient  terms. 


Using  the  constraint  equations  (3)  to  (6),  it  can  be  shown  that  four  of  the 
eleven  variables  can  be  chosen  as  solution  variables.  The  time  evolution  of  these 
solution  variables  can  be  computed  by  solving  the  corresponding  differential  equa- 
tions together  with  appropriate  initial  conditions.  The  remaining  seven  variables 
can  be  determined  from  the  solution  variables.  There  are,  however,  many  possible 
differential  equation  formulations.  The  numerical  characteristics  of  some  of  these 
formulations  will  be  discussed  in  the  next  section. 
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Table  2:  Conservative  Zone  Modeling  Differential  Equations 


Equation  Type 

Differential  Equation 

i’th  layer  mass 

pressure 

^ + qu) 

i’th  layer  energy 

i’th  layer  volume 

^ = :^((7  - 1)9.  - Kf ) 

i’th  layer  density 

dt  - CpT,V,ii^t  CplhiTt) 

i’th  layer  temperature 

dt  ~ cpPiV,(^^^  Cpth^Tt)  + ^1  dt.) 

Table  3:  Conservative  Zone  Model  Equation  Selections 


Zone  Fire  Model 

Equations 

Substitutions 

CFAST,  FAST  [5] 

d^P  dVjj  dTjj  dTj, 
dt  ^ dt  dt  I dt 

AP=P-  Pref 

CCFM.HOLE  [6] 

dAP  dy  dpu  dpi 
dt  ' dt'  dt  ' 'dt 

AP  = P - Pref,  y = VL/Aroom 

CCFM.VENTS  [4] 

dAP  dy  dmjj  dmr. 
dt  ' dt'  dt  ' dt 

AP  = P - Pref,  y = VllAroom 

FIRST,  HARVARD  V [7] 

dEjj  dEj.  dmrr  dmr. 
dt  ' dt  ' dt  ' dt 

3.2  Numerical  Characteristics  of  Several  Zone  Fire  Modeling  Dif- 
ferential Equations 

There  are  330  different  ways  to  select  four  variables  from  eleven  to  form  a system  of 
differential  equations.  Many  of  these  systems  are  incomplete  due  to  the  relationships 
that  exist  between  the  variables  given  in  equations  (3)  to  (6).  For  example  the 
variables,  pjj,  Vu,  mu  and  P form  a dependent  set  since  pu  = mufVu-  Table  3 
shows  the  solution  variable  selection  made  by  a few  zone  fire  models.  The  variable 
y that  appears  in  this  table  is  the  height  of  the  upper  layer  above  the  compartment 
floor. 

The  number  of  differential  equation  formulations  can  be  considerably  reduced 
by  not  mixing  variable  types  between  layers;  that  is,  if  upper  layer  mass  is  chosen 
as  a solution  variable,  then  lower  layer  mass  must  also  be  chosen.  For  example,  for 
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two  of  the  solution  variables  choose  mi,  and  mu,  or  pi,  and  pu,  or  Tl  and  Tf/.  For 
the  other  two  solution  variables  pick  El  and  Eu  or  P and  Vl  or  P and  Vu-  This 
reduces  the  number  of  distinct  formulations  to  nine.  Since  the  numerical  properties 
of  the  upper  layer  volume  equation  are  the  same  as  a lower  layer  one,  the  number 
of  distinct  formulations  can  be  reduced  to  six. 

The  next  several  subsections  discuss  the  numerical  implications  of  using  these 
formulations.  Some  of  the  problems  discussed  can  be  solved  by  ignoring  the  pressure 
equation.  The  resulting  approximate  equations  and  their  implications  are  discussed 
in  section  3.3. 

3.2.1  Pressure 

Some  of  the  numerical  problems  that  arise  in  zone  fire  modeling  are  due  to  the 
difficulty  of  computing  accurate  pressure  differences  across  vent  openings.  When 
adjacent  room  pressures  are  close,  a catastrophic  cancellation  will  lead  to  a loss 
of  significant  digits.  If  the  pressures  are  too  close,  the  result  of  the  subtraction  is 
roundoff  noise.  This  can  cause  problems  if  the  noise  is  amplified  in  the  next  stage 
of  the  computation;  that  is,  the  noise  is  propagated  and  may  dominate  some  later 
stage  of  the  calculation.  This  problem  is  compounded  by  the  fact  that  the  base 
pressure  in  a room  at  1 atmosphere  is  about  10^  pascals  (Pa).  Pressure  drops  across 
a vent  as  low  as  .1  Pa  can  cause  significant  mass  flow  through  a vent.  Therefore, 
if  pressure  is  used  as  a solution  variable,  seven  accurate  digits  must  be  carried  in 
order  to  have  one  significant  digit  in  the  vent  flow  calculation.  One  way  around  this 
problem  is  to  solve  for  an  offset  pressure,  AP,  where  for  some  reference  pressure, 
Pre/,  the  room  pressure  is  given  by 

P = PreJ  + AP  . 

The  differential  equation  for  relative  pressure,  AP  is  given  by 

dAP 
dt 


3.2.2  Internal  Energy 

The  problem  with  using  internal  energy  in  a formulation  is  the  difficulty  in  accurately 
determining  the  offset  pressure,  AP,  needed  for  accurate  vent  flow  calculations.  In 
large  part  this  is  due  to  the  fact  that  energy  is  not  an  intensive  variable,  it  is 
proportional  to  layer  volume  as  well  as  temperature.  To  illustrate,  consider  that 
the  total  room  pressure  can  be  expressed  in  terms  of  lower  and  upper  layer  internal 
energy  and  room  volume  by  using  equations  (3)  to  (6)  to  obtain 

P = ^{El  + Ei,). 


dP 


ref 


dt  dt 

dP 

-r-  if  Pref  is  Constant 
dt  ■’ 
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Substituting  El  = Ere/  + ^El,  Eu  = Eref  + ^Eu  and  P = Pref  + AP  into  the 
above  equation  and  solving  for  Eref  and  AP  we  obtain 


AP 


Eret 


l-^(AEi.  + AEu) 
P,dV 
2(7-1)  ■ 


The  term,  AP  can  be  small  while  the  term  \AEl\  + |AP(/|  is  large.  This  will  result 
in  cancellation  errors.  The  term,  AP,  will  not  be  zero  in  general,  therefore  we  can 
not  assume  that  AEl  = —AEu.  Section  3.3.2  discusses  the  approximation  where 
AP  is  assumed  to  be  zero. 


3.2.3  Temperature  and  Density 

The  temperature  and  density  differential  equations  have  several  advantages  and 
disadvantages  in  common.  As  seen  in  Table  2,  both  the  density  and  temperature 
equations  have  the  term  qi  — CpihiTi.  The  vent  flow  component  of  this  term  is 
identically  zero  for  flows  leaving  a zone  since  the  enthalpy  flow  rate  for  a vent 
flow  is  9vent  = CpThyent^vent-  An  Unnecessary  subtraction  can  be  avoided  by  setting 
this  vent  flow  component  to  zero  analytically,  thereby  avoiding  a loss  of  significant 
digits.  This  property  of  the  density  and  temperature  equation  is  due  to  the  fact 
that,  ignoring  pressure  transients,  the  temperature  of  a room  is  not  affected  by  flows 
leaving  it.  Similar  cancellations  must  also  be  eliminated  in  the  pressure  equation 
for  this  strategy  to  be  useful.  Unfortunately,  this  is  not  possible  unless  the  pressure 
equation  is  dropped  from  the  equation  set. 

The  density  and  temperature  equations  also  have  a problem  in  common.  Both 
have  layer  volume  terms  in  the  denominator  which  may  vanish.  Although  these 
singularities  are  removable  (if  they  were  not  the  derivatives  ^ and  ^ would  be 
infinite),  they  cause  numerical  problems  since  it  is  difficult  to  determine  proper 
values  for  ^ and  ^ in  this  case.  One  method  used  to  solve  this  problem  is  not  to 
allow  layers  to  vanish. 

3.2.4  Mass 

The  mass  equation  does  not  have  the  vanishing  denominator  problem  of  the  density 
or  temperature  equation.  Using  the  mass  equation  allows  sensible  initial  conditions. 
The  mass  for  a layer  with  zero  height  is  just  zero.  However,  quantities  derived  from 
mass  such  as  density  are  only  valid  when  a layer  volume  has  a significant  number 
of  accurate  digits. 

3.3  Approximate  Formulations 

The  formulations  discussed  in  this  section  are  approximate  in  the  sense  that  cer- 
tain terms  deemed  negligible  are  removed  from  the  modeling  differential  equations 
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Table  4:  Approximate  Zone  Model  Equation  Selections 


Zone  Fire  Model 

Equations 

Substitutions 

ASET  [8,  9] 

y — ^l!  ^Toom 

BRI  [10,  11] 

o 

II 

LAVENT  [12,  13] 

P — P . ^ 'T^  — 'T  , 

eimb?  df  ^ (If  ^ L — amo 

y — ^l!  -^room 

derived  in  Section  3.1.  Three  types  of  approximations  and  their  error  behavior  are 
discussed.  Two  involve  the  elimination  of  the  pressure  transient,  (^),  terms.  A 
third  way  approximates  the  differential  equations  by  assuming  that  the  conditions 
in  the  lower  layer  remain  at  ambient.  Some  fire  models  along  with  their  variable 
choices  that  use  some  of  the  approximation  techniques  discussed  in  this  section  are 
listed  in  Table  4.  ASET  assumes  the  pressure  remains  at  ambient  and  uses  a non- 
dimensional  form  of  the  layer  height  and  upper  layer  density  equations.  BRI  assumes 
that  the  pressure  relaxes  to  quasi-steady  values  instantly.  BRI  solves  a non-linear 
algebraic  equation  equivalent  to  equation  (9)  with  ^ set  to  zero.  First,  we  examine 
the  behavior  of  pressure  as  motivation  for  making  these  approximations. 


3.3.1  Behavior  of  the  Pressure  Equation  for  a Heated  Enclosure  with  a 
Small  Leak 


The  pressure  in  a compartment  approaches  steady  state  rapidly  if  other  compart- 
ment properties  such  as  layer  temperatures,  fire  size,  and  vent  sizes  are  constant. 
The  equilibrium  pressure  value  depends  on  the  fire  size  and  vent  areas  or  more 
generally  on  the  sources  and  sinks  of  enthalpy  in  the  room.  To  characterize  the 
equilibrium  pressure  value  and  the  time  required  to  reach  equilibrium,  consider  a 
room  with  a fire  that  is  vented  to  the  outside  as  illustrated  in  Figure  3. 

To  simplify  the  analysis  assume  that  the  vent  is  a slit  located  at  the  floor  so  that 
Bernoulli’s  law. 


^vent  — C vent  A vent  vent  ^ P ? 


can  be  used  to  model  the  mass  flow  through  the  vent.  The  conclusions  found  here 
hold  for  more  general  vent  algorithms  since  they  aU  use  or  compute  pressure  dif- 
ferences. Here  Ayent  denotes  the  area  of  the  vent,  Cyent?  the  vent  coefficient,  while 
Pvent  denotes  the  density  of  the  vent  flow  gas.  Further,  assume  that  the  density  and 
temperature  of  the  gas  flowing  through  the  vent  is  constant  over  the  time  period 
required  to  reach  pressure  equilibrium.  This  is  reasonable  since  this  time  period  is 
typically  rather  short. 
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Figure  3:  One  Room  Test  Case 


The  source  of  enthalpy  in  the  room  is  the  fire  and  is  denoted  by  9fire^*  The 
enthalpy  flow  rate  out  of  the  room  is  due  to  the  vent  and  is  denoted  by  ^vent-  The 
initial  value  problem  for  the  pressure  drop  across  the  vent  is 


dAP 

dt 

AP(0) 


■y  ( 9fire  9vent ) 

0 


(13) 


where  AP  is  the  pressure  drop  across  the  vent,  V is  the  volume  of  the  room,  and 
the  other  terms  have  been  defined  previously.  Setting 


a = 

b = 

problem  (13)  simplifies  to 


7 - 1 . 

Y 9fire 
7-1 


V 


CpCvent  •'4- ventTvent 


vent  ? 


dAP 

dt 

AP(0) 


a-b\fAP 

0 . 


This  differential  equation  is  separable  and  can  be  integrated  to  obtain  the  implicit 
solution 

‘ ^ (i-  VaV/apJ  - 

^Most  entrainment  models  used  by  fire  plume  models  do  not  affect  the  calculation  of  the  pressure 
rise  or  transient  time  since  the  enthalpy  entrained  from  the  lower  layer  cancels  with  the  enthcdpy 
added  back  into  the  upper  layer. 
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Figure  4:  Pressure  Equilibrium  Values  for  a One  Room  Case 


or  equivalently 

AP  = APoo  (15) 

with  the  equilibrium  pressure  APo©  given  by 

ap.=  (2)^i.ooxio-(^)\  (16) 

where  Tvent  = 300ii',/>vent  = l-2fcfif/m^,Cvent  = -68  and  Cp  = lOOOkJ/kgK . Taking 
the  limit  as  t — >•  oo,  equation  (15)  verifies  that  AP©©  is  the  equilibrium  pressure  and 
equation  (14)  shows  how  long  it  takes  to  achieve  it.  Substituting  AP/AP©©  = .99 
into  equation  (14)  gives  the  time  required  for  the  pressure  to  reach  99  per  cent  of 
APoo  or 


<.99  « 


21.5 


AP©© 

9fire/k^ 


For  a room  with  volume  18m^,  a 100  kW  fire,  and  an  equilibrium  pressure  of  .1 
Pa,  the  time  required  to  reach  99  per  cent  of  the  equilibrium  value  is  approximately 
.00038  seconds. 

Figure  4 shows  the  dependence  of  equilibrium  pressure  on  fire  size  and  vent  area. 
Figure  5 shows  how  the  time  required  to  reach  equilibrium  is  affected  by  these  same 
two  quantities. 

This  analysis  shows  that  pressures  tend  towards  equilibrium  values  very  quickly. 
The  equilibrium  pressure  value  is  a function  of  quantities  that  affect  the  sources 
or  sinks  of  enthalpy  in  a room.  The  flow  of  enthalpy  through  a vent  into  and  out 
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Figure  5:  Pressure  Equilibrium  Times  for  a One  Room  Case 


of  a room  is  affected  by  vent  area,  layer  heights  and  temperatures  etc.  The  fire 
also  contributes  to  the  enthalpy  gain  in  a room.  As  these  quantities  (vent  area, 
layer  heights,  fire  size,  etc)  change,  the  room  pressure  adjusts  almost  instantly  to  a 
new  quasi-steady  state  value.  Numerically,  this  property  is  known  as  stiffness.  A 
numerical  challenge  of  zone  fire  modeling  is  to  determine  pressures  accurately  and 
efficiently. 

For  a one  room  model  the  algebraic  equation  (16)  can  be  used  to  determine  the 
equilibrium  pressure  instead  of  the  pressure  differential  equation  (13).  Furthermore, 
the  ^ term  in  the  other  zone  modeling  differential  equations  (as  expressed  in  Table 
2)  can  be  dropped  in  cases  when  ^ goes  to  zero  almost  instantly. 

3.3.2  Constant  Pressure  Approximation 

The  constant  pressure  approximation  is  based  on  the  assumption  that  ^ = 0 and 
that  the  compartment  pressure  is  given  by  the  pressure  at  some  particular  elevation 
usually  the  floor.  It  is  further  assumed  that  Pqoot  is  hydrostatically  related  to  a 
reference  pressure  using 


■Pfloor  — -fref  PrefdV&ooT 

where  Pref  a.nd  pref  are  the  reference  pressure  and  density,  g is  the  acceleration 
of  gravity,  and  y^oor  is  the  distance  between  the  floor  and  the  reference  elevation. 
Usually,  Pref  is  the  atmospheric  pressure  at  the  reference  elevation  under  ambient 
conditions  and  pref  is  the  ambient  density.  In  the  simplest  case,  when  each  room’s 
floor  is  at  the  reference  elevation,  the  floor  pressure  is  the  same  as  the  reference 
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Table  5:  Approximate  Zone  Modeling  Differential  Equations 


Equation  Type 

Differential  Equation 

i’th  layer  energy 

dE,  ^ 

dt  7 

i’th  layer  volume 

rfV.  _ (7-1)9. 
dt  -yP 

i’th  layer  density 

il 

1 

1 

i’th  layer  temperature 

dTi  _ q,-Cpm,Ti 

dt  ~ Cpp,V, 

pressure  for  aU  rooms.  Table  5 gives  the  approximate  zone  modeling  differential 
equations  corresponding  to  the  conservative  ones  given  previously  in  Table  2.  The 
equations  in  this  table  were  obtained  by  neglecting  the  dP/dt  term  in  Table  2 where 
P is  the  total  pressure  in  a room. 

These  approximate  equations  are  reasonable  if  the  true  pressure  offsets  (with 
respect  to  Pfioor)  small  compared  to  the  hydrostatic  head,  the  pressure  drop 
between  the  floor  and  the  ceiling.  This  pressure  drop  is  about  35  Pa  for  a room 
with  a 3m  ceiling  containing  air  at  ambient  conditions  (p  = 1.2kg /m^).  Equation 
(16)  shows  that  for  a simple  one  room  model  an  equilibrium  pressure  offset  of  .1  pa 
corresponds  to  a ^fu-e  of  100  kw,  an  Avent  of  Im^,  and  ambient  pvent  and  Pvent- 

The  absolute  errors  generated  by  using  the  approximate  equations  in  Table  5 
can  be  estimated  by  defining  an  error  differential  equation  for  each  variable.  As- 
suming that  the  mass  and  enthalpy  sources  are  the  same  for  the  conservative  and 
approximate  differential  equations,  define  the  approximation  error  for  the  tempera- 
ture equation  using 

where  T,  denotes  the  approximate  temperature.  An  initial  value  problem  for  ef  can 
be  derived  using 


de^ 

t 

dfi 

dTi 

dt 

dt 

dt 

-1 

dP 

CpPi 

dt 

ef(0)  = 

0 . 

This  equation  can  be  solved  assuming  that  the  density  is  constant  over  the  time 
period  required  to  achieve  equilibrium.  The  result  is 

-APc 

temperature  error  ■ 

CpPi 
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Table  6:  Approximate  Zone  Modeling  Differential  Equation  Error 


Equation  Type 

Error  Estimate 

i’th  layer  energy 

7 

i’th  layer  volume 

VijAPpo 

7P 

i’th  layer  density 

APpn 

CpTi 

i’th  layer  temperature 

-APoo 

Cppi 

The  error  in  temperature  resulting  from  this  approximation  is  related  to  the  mag- 
nitude of  the  equilibrium  pressure.  The  approximate  equations  are  then  valid  eis 
long  as  the  pressure  rise  is  not  significant.  Another  way  of  looking  at  this  is  that 
the  conservative  equations  in  Table  2 are  essentially  the  same  as  the  approximate 
equations  in  Table  5 except  during  a transient  period  when  ^ decays  to  zero. 

Approximate  errors  for  the  other  zone  modeling  equations  are  derived  similarly 
and  are  presented  in  Table  6.  The  error  defined  by  equation  (17)  underestimates 
the  actual  error  since  the  enthalpy  flow  rate  out  of  a room  is  highly  sensitive  to  the 
pressure  drop  across  the  vent. 


3.3.3  Algebraic  Pressure  Approximation 

The  pressure  differential  equation  (9)  can  be  approximated  algebraically  by  deter- 
mining the  pressure  that  satisfies  the  non-linear  equation 

iL(^)  + = 0 . (18) 


Unlike  the  pressure  approximation  in  the  previous  section,  this  pressure  will  not  in 
general  be  constant  in  time.  For  simple  cases,  such  as  a one  room  model  with  one 
fire  and  layer  height  above  a thin  slit  vent,  this  equation  can  be  solved  analytically 
for  F.  For  example,  using  equation  (16)  for  the  equilibrium  pressure  offset  AFqo, 
this  equation  has  the  solution 


^ — ^ref  + 


Qfire 


CpCyent  A vent^vent  \/^Pvent 


2 


Rehm  and  Baum  in  [14]  similarly  calculate  the  equilibrium  pressure  rise  in  the 
context  of  field  modeling. 

For  multiple  room  simulations  the  approximation  is  more  complicated.  An  equi- 
librium pressure  offset  must  be  calculated  for  each  room  using  an  equation  similar 
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to  (18).  This  results  in  a system  of  equations  since  the  equilibrium  pressure  in  one 
room  depends  on  the  equilibrium  pressure  in  another  via  vent  connections.  Differen- 
tial equation  solvers  such  as  Petzold’s  DASSL  [15,  16]  can  solve  algebraic  equations 
simultaneously  with  differential  equations.  This  was  tried  with  CCFM. VENTS. 
It  was  found  that  there  was  no  advantage  to  doing  this.  The  algebraic  pressure 
equation  version  could  not  track  rapidly  changing  pressures  and  run  times  were  not 
significantly  shorter. 

3.3.4  A One  Zone  Approximation 

Conditions  in  the  lower  layer  are  essentially  ambient  for  many  fire  scenarios  of 
interest.  This  observation  may  be  exploited  by  replacing  the  differential  equation 
for  lower  layer  density  (equation  (11))  with  pi,  = Pamb  or  the  differential  equation 
for  lower  layer  temperature  (equation  (12))  with  Ti  = Tamb-  This  is  equivalent  to 
assuming  that  the  rate  of  mass  and  enthalpy  additions  to  the  lower  layer  satisfy 
qL/{cpmL)  = Tamb  which  wiU  be  true  as  long  as  the  temperature  of  flows  added  to 
the  lower  layer  are  at  ambient. 

4 Zone  Fire  Modeling  Numerical  Characteristics 

Physical  models  of  natural  and  forced  vent  flow,  fire  plumes,  radiation,  conduction, 
convection  etc  are  used  to  exchange  or  transfer  mass  and/or  energy  between  zones. 
As  with  the  differential  equation  formulations  discussed  in  the  previous  section,  two 
physical  models  for  computing  these  phenomena  may  be  identical  physically  but  be 
different  numerically.  This  section  addresses  some  of  the  numerical  issues  important 
in  the  design  of  physical  algorithms,  some  miscellaneous  numerical  considerations 
appropriate  for  any  fire  model  and  finally  discusses  the  numerical  properties  of 
systems  of  differential  equations. 

4.1  Numerical  Characteristics  of  Some  Physical  Models 

4.1.1  Natural  Vent  Flow 

Numerical  difficulties  can  arise  when  calculating  vent  flow  because  of  its  dependence 
on  the  square  root  of  pressure  differences.  This  is  especially  a problem  when  the 
pressure  drop  across  the  vent  is  small  relative  to  the  pressure  computed  in  each  room 
adjacent  to  the  vent.  To  illustrate  this  phenomena  consider  two  rooms.  Suppose 
the  first  room  has  a fire  and  is  connected  to  the  outside  and  to  a second  room  which 
does  not  have  a fire.  This  second  room  is  assumed  to  be  connected  only  to  the  first 
room  and  not  to  the  outside.  To  simplify  the  analysis  assume  that  the  vents  are 
narrow  slits  located  at  the  floor.  This  configuration  is  depicted  in  Figure  6. 

Suppose  that  the  pressure  above  ambient  in  rooms  1 and  2 are  denoted  APi 
and  AP2.  Theoretically,  these  pressure  offsets  will  be  the  same  after  the  initial 
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Figure  6:  Two  Room  Test  Case  Configuration 


pressure  transient  dies  away  according  to  the  analysis  in  Section  3.3.1.  Numerically, 
however,  these  two  offsets  will  be  different.  Though  this  difference  may  be  insignif- 
icant physically,  the  magnification  of  the  cancellation  error  occurring  in  the  vent 
flow  computation  may  be  significant  numerically.  Unfortunately,  simply  setting  the 
pressure  drop  to  zero  for  physically  insignificant  flows  is  not  an  adequate  solution 
to  this  problem  of  unwanted  error  propagation  as  will  be  explained  next. 

In  this  test  problem,  the  pressure  in  room  1 will  rise  to  its  equilibrium  value 
as  given  by  equation  (16)  in  a time  period  given  by  equation  (14).  If  a cut  off 
pressure  drop  APcut  is  chosen  so  that  = 0 for  IAP2  ~ APi|  < APcut?  then  APi 
will  continue  to  rise  (since  is  not  zero)  while  AP2  remains  fixed.  Eventually 
IAP2  — APil  win  exceed  APcut  so  that  will  not  be  zero.  Again,  AP2  wiU  rapidly 
approach  APi  until  IAP2  — APi|  < APcut-  This  type  of  algorithm  will  result  in  a 
drastic  reduction  in  time  step  size  as  the  differential  equation  solver  tries  to  track 
the  solution  due  to  the  “stair  stepping”  behavior  of  AP2  illustrated  in  Figure  7.  Of 
course,  for  this  test  problem  special  methods  may  be  used  to  resolve  the  problem 
such  as  eliminating  the  differential  equation  for  by  setting  AP2  = APi.  This 
will  not  work  for  the  general  case. 

Denote  the  pressure  difference  across  the  vent  by  AP  = APj  — AP2.  This  section 
establishes  a criterion  for  choosing  APcut  such  that  pressure  differences  satisfying 
AP  > APcut  contain  at  least  some  accurate  digits  and  gives  a method  for  smoothly 
damping  AP  to  zero  when  AP  < APcut* 

Although  the  Bernoulli  law  for  computing  vent  flow  is  known  to  break  down  for 
small  pressure  differences,  the  problems  discussed  here  are  solely  numerical  and  are 
a consequence  of  the  fact  that  only  a finite  number  of  digits  are  used  to  compute 
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Figure  7:  Numerical  Pressure  Behavior  for  Two  Room  Test  Case 


and  represent  the  pressure  offsets  APi  and  AP2. 

Suppose  that  A Pi  is  an  o(l)  quantity  and  is  known  to  an  absolute  error  tolerance 
of  10“^.  The  binary  representation  of  APi  then  has  about  10  base  2 digits  of 
information.  This  number  is  obtained  by  solving  the  equation  2“”  = 10“^  for  n. 
The  other  14  base  2 digits  (assuming  a 24  bit  mantissa)  contain  roundoff  noise. 
If  the  first  13  base  2 digits  of  two  o(l)  pressure  offsets  APi,  AP2  agree,  then  the 
subtraction  APi  — AP2  will  result  in  loss  of  aU  significant  digits.  It  is  therefore 
important  to  know  when  a subtraction  will  result  in  total  cancellation. 

The  number  of  decimal  digits  lost  in  the  subtraction  APi  — AP2  can  be  estimated 
using 

majc(l,  APi,  AP2) 

Wlost  - ogio  |APi  - AP2I 

If  logio  is  replaced  by  log2  then  niost  estimates  the  number  of  binary  digits  lost. 
Suppose  that  APi  and  AP2  each  have  m significant  digits  and  suppose  that 


niost  < rn  . 

Substituting  equation  (19)  and  exponentiating,  it  follows  that 


lAPi  - AP2I  > Cpmax(l,  APi,  AP2) 


where  Cp  = 10“”^.  This  inequality  gives  a criterion  for  deciding  if  APi  - AP2  has  any 
accurate  digits.  To  smoothly  damp  the  roundoff  noise  present  in  the  computation 
AP  = APi  — AP2,  AP  can  be  replaced  by 


AP  = AP 
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where  Ccut  defined  by 


Ccut  = 10€pmax(l,  APi,  AP2)  • 

When  AP  is  large  relative  to  Ccut?  a-nd  AP  are  essentially  the  same.  On  the 
other  hand,  when  AP  is  small  relative  to  tcut?  AP  is  essentially  zero. 

4.1.2  Forced  Vent  Flow 

A forced  ventilation  system  consists  of  a network  of  ducts  and  supply  and  exhaust 
fans.  In  general,  the  volume  flow  rate  delivered  by  a fan  depends  on  the  pressure 
drop  across  the  fan.  This  relation  is  known  as  the  fan  curve.  Numerical  problems 
can  occur  even  in  simple  systems  consisting  of  one  fan  with  a constant  volume  flow 
rate  and  one  duct  connecting  two  rooms.  For  example,  consider  a fan  exhausting 
a room  where  sources  of  enthalpy  such  as  fires  and  natural  vents  are  small  relative 
to  the  enthalpy  removed  by  the  exhaust  fan.  Numerical  problems  occur  when  the 
fan  is  modeled  by  removing  flow  at  a single  elevation.  In  this  case,  the  fan  is  either 
in  the  lower  layer  or  the  upper  layer.  While  the  fan  is  in  the  lower  layer,  the  layer 
interface  drops  since  the  fan  is  “stronger”  than  the  fire.  Similarly,  while  the  fan 
is  in  the  upper  layer  the  layer  interface  rises.  This  problem  is  ill-conditioned  since 
the  total  enthalpy  flow  rate,  {qi  + qu),  is  a discontinuous  function  of  layer  height. 
When  the  layer  interface  is  near  the  fan  elevation,  small  changes  in  the  layer  height 
produce  large  changes  in  the  layer  height  derivative.  A differential  equation  solver 
detecting  this  jump  in  the  layer  height  derivative  (or  layer  volume  derivative)  will 
reduce  the  time  stepsize  drastically  to  try  and  track  the  layer  height. 

The  solution  to  this  problem  is  to  model  the  fan  by  withdrawing  air  over  a finite 
but  nonzero  height  in  the  wall.  When  the  layer  is  within  a forced  fan  duct  flow  is 
withdrawn  from  both  layers.  This  results  in  smoother  transition  in  the  layer  height 
derivative. 


4.1.3  Fire  Plumes 


A fire  deposits  energy  denoted  q^^,  q^^  into  the  lower  and/or  upper  layer.  This 
lower  layer  contribution  is  set  to  zero  in  most  plume  models.  A fire  also  transfers 
mass  and  energy  from  the  lower  layer  to  the  upper  layer  by  entraining  a portion 
of  the  lower  layer  into  the  hot  plume  which  in  turn  flows  due  to  buoyancy  into  the 
upper  layer.  Let  the  energy  entrainment  terms  be  denoted  and  The  total 
energy  contribution  to  the  lower  and  upper  layers  due  to  both  mechanisms  is  then 

Qplume  = 9L  + 9^nt.  (20) 

9^ume  = + (21) 

It  is  usually  assumed  that  all  of  the  energy  taken  out  of  the  lower  layer  due  to 
entrainment  is  added  back  into  the  upper  layer  so  that  g^^  = — g^f  Therefore 


•L  -u  - aL  , aU 

Vplume  ' Tfplume  9fire  ' 9fire 


(22) 
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To  avoid  a catastrophic  cancellation,  a plume  model  should  calculate  the  total  energy 
contribution  using  equation  (22)  rather  than  summing  the  terms  as  calculated  in 
equations  (20)  and  (21).  This  is  especially  important  when  is  large  compared 
9fire*  total  energy  term  is  required  by  zone  fire  models  which  use  the 
pressure  differential  equation  or  its  equivalent. 

4.2  Numerical  Characteristics  of  the  Zone  Fire  Modeling  Differen- 
tial Equations 

In  Section  3.2  the  numerical  properties  of  individual  differential  equations  were 
discussed.  It  is  important,  however,  to  understand  how  these  equations  behave 
together  as  a system  in  order  to  implement  accurate  and  efficient  algorithms  for 
their  solution.  The  key  property  that  zone  modeling  differential  equations  possess 
is  called  stiffness.  Differential  equation  solvers  not  taking  this  property  into  account 
wiU  at  best  be  grossly  inefficient  and  at  worst  give  wrong  answers. 

4.2.1  Comparison  of  the  Physical  and  Numerical  Behavior  of  the  Dif- 
ferential Equations 

Stiff  systems  of  ordinary  differential  equations  (ODE’s)  are  an  important  class  of 
problems  that  can  occur  when  the  modeled  phenomenon  possesses  characteristic 
time  scales  that  vary  by  several  orders  of  magnitude.  Physically,  the  system  of 
ODE’s  used  in  zone  fire  modeling  are  stiff  because  the  pressure  adjusts  to  changing 
conditions  much  faster  than  other  quantities  such  as  upper  and  lower  layer  temper- 
atures or  layer  heights. 

The  numerical  difficulties  encountered  because  of  stiffness  can  not  be  avoided 
by  exchanging  the  pressure  equation  for  some  other  equation  such  as  temperature, 
density,  or  internal  energy.  As  shown  in  Table  2,  each  zone  modeling  differential 
equation  contains  a ^ term.  If  the  pressure  is  computed  using  one  of  the  approx- 
imations discussed  in  sections  3.3.2  or  3.3.3  and  ^ is  removed  from  the  modeling 
differential  equations,  the  resulting  approximate  ODE’s  are  not  stiff  and  a standard 
nonstiff  solver  may  be  used.  However,  the  class  of  problems  that  can  be  solved  is 
reduced  since  large  pressure  fluctuations  can  not  be  modeled  properly. 

The  curious  aspect  of  stiff  ODE’s  is  that  the  solution  appears  to  be  chang- 
ing slowly  and  yet  the  computational  costs  of  computing  the  solution  are  enormous 
when  using  standard  nonstiff  algorithms  such  as  Runge-Kutta  methods  or  predictor- 
corrector  methods  such  as  Adams-Bashforth.  The  question  then  is  why  does  it  cost 
so  much  to  solve  a problem  whose  solution  changes  slowly?  To  maintain  stability, 
a nonstiff  solver  must  use  a stepsize  that  is  small  enough  to  track  the  part  of  the 
solution  corresponding  to  the  shortest  time  scale  even  when  this  solution  component 
decays  rapidly  to  some  quasi-steady  value.  This  stepsize  is  much  smaller  than  re- 
quired to  accurately  track  the  desired  part  of  the  solution  which  corresponds  to  one 
of  the  longer  time  scales.  So  for  stiff  problems  the  choice  of  stepsize  is  dominated 
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by  considerations  of  stability,  not  accuracy.  These  concepts  and  their  significance 
will  be  discussed  in  section  4.2.2. 

Stiff  differential  equation  solvers  are  expensive  to  use.  This  is  because  a nonlinear 
set  of  simultaneous  equations  involving  the  solution  variables  must  be  solved  at  each 
time  step.  Newton’s  method  is  typically  used  for  this  and  requires  the  solution  of  a 
system  of  linear  equations  at  each  Newton  iteration.  Consequently,  it  is  inefficient 
to  use  stiff  methods  for  nonstiff  problems.  Stiff  methods  work  because  their  choice 
of  stepsize  is  dominated  by  considerations  of  accuracy,  not  stability.  As  a result, 
they  allow  larger  time  steps  to  be  taken  than  nonstiff  methods.  Even  though  the 
work  per  step  is  greater,  the  number  of  steps  is  sufficiently  smaller  that  the  total 
work  is  smaller. 

Several  approaches  can  be  taken  to  decide  whether  a given  system  of  differential 
equations  is  stiff.  One  approach  is  to  let  the  solver  decide.  Many  modern  ODE 
solvers  have  built-in  heuristics  to  determine  if  a problem  is  stiff.  To  analyze  the  nu- 
merical character  of  the  ODE’s  in  the  zone  fire  model  CCFM. VENTS,  DEPAC  was 
used.  DEPAC  is  a set  of  three  ODE  solvers  DERKF,  DEABM,  and  DEBDF  designed 
by  Shampine  and  Watts.  DERKF  is  a fifth  order  variable  step  size  Runge-Kutta 
code.  It  can  be  used  for  nonstiff  and  mildly  stiff  ODE’s  when  derivative  evaluations 
are  not  expensive.  It  should  not  be  used  for  high  accuracy,  nor  for  answers  at  a great 
many  points.  DEABM  is  a variable  order,  variable  step  size  Adams  code.  It  can  be 
used  for  nonstiff  and  mildly  stiff  ODE’s  when  high  accuracy  is  required.  DEBDF  is 
a variable  order,  variable  step  size  backward  differentiation  formula  code  which  can 
be  used  on  stiff  ODE’s  when  moderate  accuracy  is  required.  DERKF  and  DEABM 
attempt  to  determine  when  their  use  is  not  suitable  by  performing  diagnostics  for 
stiffness.  When  used  in  CCFM.VENTS  both  DERKF  and  DEABM  reported  that 
the  problem  might  be  stiff.  This  occurred  for  a wide  range  of  fire  scenarios. 

A second  approach  is  to  analyze  the  Jacobian  of  the  right  hand  side  of  the  sys- 
tem of  ODE’s.  A subroutine,  EIGF,  was  written  to  approximate  this  Jacobian  and 
its  eigenvalues.  The  characteristic  local  time  scales  of  the  solution  were  determined 
from  the  eigenvalues.  It  was  found  that  the  characteristic  time  scales  for  the  pres- 
sure variables  were  much  smaller  than  for  the  other  solution  variables.  For  some 
problems,  the  variation  in  time  scales  was  over  five  orders  of  magnitude.  As  a result 
of  this  analysis,  the  stiff  ODE  solver  DEBDF  was  chosen. 

4.2.2  Some  Numerical  Considerations  of  Solving  Stiff  Differential  Equa- 
tions 

A zone  fire  modeling  initial  value  problem  can  be  expressed  in  vector  form  as  follows. 
At  time  t > to  find  the  value  to  the  solution  y to 

(23) 

y{to)  = yo 
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where  y and  / are  a real  valued  N-vector  functions.  The  vector  y contains  the 
solution  variables  chosen  and  the  right  hand  side  / gives  the  rate  at  which  these 
variables  change  with  time. 

Modern  ODE  solvers  use  a variable  stepsize  strategy;  that  is,  the  solver  takes  a 
number  of  steps  to  integrate  (23)  from  to  to  /,  but  these  steps  are  not  generally  all 
of  the  same  size.  Suppose  that  K steps  are  taken  and  approximations  t/,-  are  found 
to  y{T)  at  T = = t.  Consider  the  “local”  initial  value  problem:  find  the 

solution  at  t,-  to 

— it)  = /(^,^) 
z{U-\)  = Vi-i  • 

The  local  error  is  defined  by  z{ti)  — yi,  for  i = 0,.  ..,K.  Most  solvers  control  the 
local  error  at  each  step.  The  true  or  global  error  is  defined  by  e,-  = y{ti)  — yi, 
i = 0, . . . , K . Note  that  = yit)  — yK  is  the  error  in  the  approximate  solution  at 
time  t.  Consider  the  decomposition 

e,-  = {y{U)  - z{U))  + {z{ti)  - yi)  . (24) 

A linearization  shows  that  to  terms  of  first  order 

(y(r)  - z{t))'  = fy{T,  y{T)){y{T)  - z{t))  . 

The  N X N matrix  fy  is  called  the  Jacobian  matrix.  If  the  initial  value  problem 
(23)  is  stable  (all  the  eigenvalues  of  fy  have  negative  real  parts),  then  the  first  term 
in  (24)  will  not  grow  with  index  i.  The  second  term  is  the  local  error  and  will  be 
dependent  on  the  accuracy  and  stability  of  the  numerical  method  employed.  If  the 
numerical  method  is  stable,  then  the  growth  of  e,-  with  i will  be  at  most  linear  (and 
not  exponential). 

Let  yi  denote  the  numerical  solution  determined  from  exact  data  y(t,_i),  y(t,_2),. . .. 
The  following  alternative  decomposition  of  the  global  error  is  also  useful: 

local  truncation  error  propagation  error 

e.-  = (y{U)  - yi)  + (yi  - yi)  • (25) 

The  propagation  error  can  be  written  as  a product  of  an  amplification  factor  and 
the  error  term  e,_i.  If  the  numerical  method  is  stable,  the  absolute  value  of  the 
amplification  factor  will  be  less  than  or  equal  to  one  and  the  propagation  error  term 
will  not  grow  with  i.  The  local  truncation  error  depends  on  the  accuracy  of  the 
numerical  method.  These  errors  at  z = 2 are  illustrated  in  Figure  8. 

There  is  no  one  definition  of  stiffness  that  is  universally  applied  to  initial  value 
problem  (23).  One  that  is  commonly  applied  is  the  following  (see  [17]).  Definition: 
The  initial  value  problem  (23)  is  called  stiff  (oscillatory)  if  the  eigenvalues,  Xj  = 

Uj  ivj,  j = 1, . . .,  A of  the  Jacobian,  fy,  satisfy 

Uj  < 0,  i = 1,...,A  , 


26 


Solution  to  the 


Solution  to: 


dy/dt  = f(y,t) 

y(h)=  yi 


Figure  8:  Differential  Equation  Error  Terms 
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and 


mcLx(|nj|)>  niin(|nj|);. 
i<J<N  ^ l<j<N  ^ 

A physical  interpretation  of  stiffness,  which  is  imprecise  mathematically,  is  to  say 
that  the  initial  value  problem  (23)  is  stiff  if  the  process  being  modeled  possesses 
multiple  time  scales  which  vary  widely.  As  has  been  discussed  earlier,  this  is  the 
case  with  zone  fire  models  because  pressures  equilibrate  much  quicker  than  do  layer 
volumes,  masses,  densities,  or  temperatures.  Yet  another  (imprecise)  way  to  char- 
acterize stiffness  woidd  be  to  say  a problem  is  stiff  if  explicit  solvers  such  as  DERKF 
or  DEABM  cannot  handle  it.  ODE  solvers  are  either  explicit  or  implicit.  Implicit 
solvers  are  generally  stable  for  a much  wider  range  of  stepsizes  than  are  explicit 
solvers.  If  a solver  were  capable  of  performing  the  decomposition  (25),  it  would 
regard  a problem  as  stiff  if  the  stepsize  required  to  maintain  stability  (keep  the 
propagation  error  term  from  growing)  is  much  smaller  than  that  required  to  keep 
the  local  truncation  error  term  small.  In  this  case,  the  stepsize  is  said  to  be  re- 
stricted by  stability,  rather  than  by  accuracy.  Unfortunately,  there  is  no  practical 
way  to  make  this  characterization  quantitative.  What  is  estimated  at  each  step 
is  the  local  error.  Stepsize  choice  is  based  on  this  estimate.  Diagnostic  messages 
are  issued  when  the  stepsize  is  extremely  small  or  the  number  of  steps  taken  has 
exceeded  a predetermined  bound.  These  messages  typically  suggest  that  the  user 
try  a stiff  solver  such  as  DEBDF. 

Implicit  methods  are  more  suitable  for  stiff  problems.  In  the  case  of  the  backward 
differentiation  formulas  used  in  DEBDF,  the  stepsize  choice  for  a stiff  problem  will 
generally  be  restricted  by  accuracy,  not  by  stability.  Stepsize  can  be  increased  at  a 
fixed  accuracy  level  by  increasing  the  order  (of  accuracy)  of  the  method.  DEBDF 
provides  formulas  of  orders  one  through  five,  all  of  which  have  excellent  stability 
properties. 

The  simulation  time  interval  [to,  t]  for  a zone  fire  model  can  be  broken  into  two 
types  of  subintervals,  stiff  transient  and  stiff.  In  the  stiff  transient  subintervals, 
the  pressure  is  rapidly  moving  toward  a quasi-steady  state  value.  The  simulation 
generally  begins  with  a stiff  transient.  Each  time  layer  height  passes  a vent  boundary 
or  the  fire  output  takes  a jump,  a new  stiff  transient  begins.  During  a stiff  transient, 
stepsize  wiU  generally  be  small  because  it  will  be  restricted  by  accuracy;  that  is,  a 
small  stepsize  wiU  be  required  to  accurately  track  the  rapidly  changing  pressure. 
Nonstiff  solvers  can  generally  integrate  over  the  stiff  transient  time  subinterval. 
Outside  of  these  very  short  time  intervals,  a stiff  solver  is  required.  Since  there  is 
a large  overhead  associate  with  switching  solvers,  it  is  more  efficient  to  use  the  stiff 
solver  throughout  the  computation. 
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5 Conclusions 


Tools  were  presented  for  analyzing  the  numerical  behavior  of  fire  algorithms.  These 
tools  show  that  natural  vent  flow  is  particiilarly  susceptible  to  numerical  problems 
due  to  the  loss  of  significant  digits  that  can  occur  when  computing  pressure  dif- 
ferences. One  approach  for  reducing  this  problem  is  to  modify  the  subtraction 
|APi  — AP2I  in  order  to  damp  out  unwanted  error  that  occurs  when  APi  and  AP2 
are  nearly  equal. 

Many  model  formulations  can  be  derived  from  the  basic  mass  and  energy  con- 
servation equations  some  of  which  are  analytically  equivalent;  that  is,  the  equations 
in  one  formulation  can  be  converted  into  those  of  another  using  the  ideal  gas  law 
and  definitions  of  density  and  internal  energy.  These  formulations,  however,  are  not 
equivalent  numerically.  When  pressure  changes  are  significant,  the  pressure  equa- 
tion should  be  included  in  a multi-room  model.  The  significance  of  pressure  changes 
can  be  evaluated  in  terms  of  flow  through  a vent  since  vent  flow  is  quite  sensitive 
to  small  pressure  changes.  When  it  is  valid  to  assume  that  pressures  are  constant, 
the  modeling  differential  equations  can  be  simplified  by  setting  ^ to  zero.  The 
resulting  ODE’s  have  the  added  advantage  that  they  are  easier  to  solve  since  they 
are  no  longer  stiff. 
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Nomenclature 


a,  b 

•^vent 

b 

Ca(/)(Cr(/)) 

Cp(Ci,) 

Cyent 

T 

e/ 

e, 

E 

Eref 

e 

^min  (^max) 
/ 

fy 

fl() 

h 

K 

m 

J^lost 

p 

Pent 


PreJ 

P amb 

QlAqu) 

Qfiie 

Qvent 

R 


miscellaneous  constants  used  in  the  equilibrium  pressure  relation, 
equation  (16) 

area  of  vent  vn? 

base  used  in  the  Wilkinson  floating  point  model 

absolute  (relative)  condition  number  of  f 

specifle  heat  at  constant  pressure  (volume)  kJ/kg  K 

vent  coefficient,  usually  has  a value  of  around  .68 

absolute  error  of  temperature  K 

global  error  of  the  differential  equation  solution 

internal  energy  kJ 

reference  internal  energy  kJ 

Exponent  in  Wilkinson  floating  point  model 

smallest  (largest)  exponent  in  the  Wilkinson  floating  point  model 

a vector  valued  function  which  is  the  right  hand  side  of  the  dif- 
ferential equation  to  be  solved 

Jacobian  of  the  right  hand  side  vector  function  / 

The  operator  fl  maps  a real  number  to  its  floating  point  repre- 
sentation. 

Step  size 

Vent  flow  coefficient  kJ>/m/kg 
mass  kg 

number  of  digits  (decimal)  lost  due  to  cancellation  in  a subtrac- 
tion 

absolute  pressure  Pa 

pressure  below  which  flows  (through  vents)  are  deemed  to  be 
negligible  Pa 

reference  pressure  Pa 

ambient  pressure  usually  about  101325  Pa 

total  enthalpy  flow  rate  into  the  lower  (upper)  layer  W 

energy  release  rate  of  a Are  W 

energy  flow  rate  through  a vent  W 

2 

universal  ideal  gas  constant 
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T 

to 

t.99 

Tvent 

t 

t 

Tkinb 

V 

y 

y 

yi 

yo 

AP 

AP 

AP 

APi,(AP2) 

^Pcut 
A Poo 


1 

P 

Pamb 

P\ent 

T 


temperature  K 
initial  time  s 

time  for  pressure  in  a room  to  reach  99  percent  of  its  final  value 
s 

temperature  of  gas  flowing  through  vent  K 

number  of  digits  in  the  Wilkinson  floating  point  model 

independent  variable,  time,  in  the  differential  equations 

ambient  temperature 

volume 

layer  height  m 

solution  of  the  differential  equation 
numerical  solution  of  a differential  equation  at  step  i 
solution  of  a differential  equation  at  time  to 
pressure  drop  across  a vent  Pa 

pressure  drop  across  a vent  modified  to  account  for  catastrophic 
cancellation  errors  Pa 

Pressure  offset  satisfying  P = P^e/  + AP  Pa 
pressure  rise  above  floor  pressure  in  room  1 (2)  Pa 
pressure  below  which  flow  is  insignificant  numerically  Pa 
Equilibrium  pressure  Pa 

machine  epsilon,  smallest  number  added  to  1 on  the  computer 
that  gives  a result  different  than  1. 

error  tolerance  of  pressure  variables  used  in  the  calculation  vent 
flows 

ratio  of  specific  heats  Cp, 

density  kg/m^ 

ambient  density  kg/m^ 

density  of  gas  flowing  through  vent 

time  s 
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A A Description  of  Floating  Point  Arithmetic 

A model  of  floating  point  arithmetic  consists  of  a description  of  how  floating  point 
numbers  are  represented  (word  length  and  bits  for  sign,  exponent,  and  mantissa), 
how  operations  with  these  numbers  behave  (error  properties),  and  what  happens 
when  an  arithmetic  fault  occurs  such  as  division  by  zero. 

Properties  of  the  real  number  line  such  as  closure,  distributivity  and  associativity 
do  not  hold  for  “computer”  numbers.  Consequently,  two  algorithms  may  be  equiva- 
lent analytically  but  not  the  same  numerically.  The  methods  used  for  implementing 
equations  that  arise  in  Are  modeling  are  therefore  important. 

An  early  problem  in  scientiflc  computation  was  the  lack  of  a uniform  floating 
point  environment.  Each  computer  manufacturer  had  their  own  method  for  repre- 
senting and  using  floating  point  numbers.  As  a result,  scientiflc  programs  running 
on  different  computers  would  occasionally  produce  significantly  different  numerical 
results. 

Wilkinson,  considered  the  father  of  floating  point  error  analysis  describes  [18,  19] 
a simple  model  of  floating  point  arithmetic.  To  address  the  problem  of  lack  of 
uniformity  in  floating  point  implementations,  Kahan  [20,  21]  extended  the  Wilkinson 
model  and  proposed  a standard  for  floating  point  arithmetic.  Kahan’s  proposal 
along  with  features  of  several  other  proposals  became  the  IEEE  Standard  for  Binary 
Floating  Point  Arithmetic  [22,  23,  24]. 

The  IEEE  Standard  uses  three  formats  to  represent  floating  point  numbers,  sin- 
gle, double  and  extended  precision.  Single  uses  32  bits,  double  uses  64  bits.  The 
format  of  extended  arithmetic  is  left  open.  Two  special  bit  patterns  designated 
NaN  (not  a number)  and  infinity  are  used  to  represent  the  results  of  special  op- 
erations. When  the  result  of  a mathematical  operation  is  undefined  such  as  zero 
divided  by  zero  a result  of  NaN  is  returned.  Some  compilers  initialize  variables  with 
NaN  to  detect  the  use  of  undefined  variables.  The  standard  also  specifies  that  the 
operations  addition,  subtraction,  multiplication,  division,  square  root,  remainder 
and  comparison  be  provided.  The  default  rounding  behavior  takes  the  result  of 
an  operation  as  if  correct  to  infinite  precision  and  rounds  it  to  the  nearest  machine 
number.  Rounding  the  result  toward  or  away  from  zero  is  also  provided.  Rounding 
towards  zero  is  sometimes  called  chopped  or  truncated  arithmetic. 

The  real  number  line  is  modeled  on  the  computer  as  a finite  collection  of  floating 
point  numbers  where  each  floating  point  number  consists  of  a mantissa  or  fraction, 
an  exponent,  and  a sign.  These  numbers  are  sometimes  called  model  or  machine 
numbers.  Wilkinson  [18,  19]  parameterizes  the  mantissa  and  exponent  in  terms  of 
the  four  parameters:  6,  t,  Cmin  and  Cmax  where  b represents  the  base  (6  = 2 for  binary, 
6 = 16  for  hexa-decimal),  t is  the  precision  or  number  of  digits  in  the  mantissa,  Cnun 
is  the  smallest  exponent  and  e^ax  is  the  largest.  Values  of  these  parameters  for 
various  computers  are  presented  in  Table  7. 

The  IEEE  Standard  is  a model  of  binary  floating  point  arithmetic.  Most  comput- 
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Table  7:  Single  Precision  Floating  Point  Characteristics  for  Various  Computers 


Computer 

b 

t 

€min 

^max 

€ - b^-^ 

Uni  vac  1100 

2 

27 

-128 

127 

1.49  X 10-* 

PDP-ll/Vax  11/780 

2 

24 

-127 

127 

1.19  X 10“'^ 

Cray  - 1 

2 

47 

-8189 

8190 

1.42  X lO-i'* 

Apple  Macintosh,  IBM  PC 

2 

24 

-128 

127 

1.19  X lO-"^ 

IBM  360/370,  Concurrent  7/32 

16 

6 

-64 

63 

9.54  X lO-"^ 

CDC  200  (Cyber  205) 

2 

47 

-28625 

28718 

1.42  X 10-1^ 

IEEE  754  Standard 

2 

24 

-128 

127 

1.19  X 10"^ 

ers  built  today  use  this  standard  or  make  it  available  as  an  option.  The  Apple  Mac- 
intosh, the  IBM  PC  and  clones,  and  UNIX  workstations  from  Sun,  Silicon  Graphics, 
DEC,  IBM,  and  HP  aU  use  the  IEEE  Standard  for  floating  point  arithmetic.  Im- 
plementations differ  from  vendor  to  vendor,  but  these  vendors  claim  to  adhere  to 
the  IEEE  Standard.  Other  manufacturers  such  cis  Convex  make  it  available  as  an 
option.  The  Cray  series  of  computers  is  a notable  exception.  Programs  that  run 
on  computers  that  use  IEEE  Standard  arithmetic  should  give  similar  floating  point 
results. 

The  floating  point  numbers  of  the  Wilkinson  model,  on  which  the  IEEE  Standard 


is  based,  can  be  generated  from  the  four  parameters  6,  t,  Cnum 
following  representation 

Cmax  using  the 

± (6i  X 6“  + 62  X 6-'  + ...  + 6,  X 6‘-‘)  X 6' 

(26) 

where 

1 

VI 

VI 

(27) 

0 < bi  < b — 1,  i = 2, . . .,t 

(28) 

^min  ^ C ^ ^max* 

(29) 

The  numbers  6i,. . bt  form  the  mantissa  and  e denotes  the  exponent.  If  the 
first  bit  of  the  mantissa,  6i,  is  never  0 then  the  set  of  floating  point  numbers  is 
said  to  be  normalized.  For  normalized  binary  systems  the  first  bit  is  usually  not 
stored  (since  it  is  always  1).  Then  32  bits  can  then  be  used  to  represent  1 sign  bit, 
an  8 bit  exponent  and  a 24  bit  mantissa. 
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b=2 

t=3 

emm=-l 
emax  = 2 


c 

1/ 

2 

III 

1 

1 1 1 

2 

1 1 

4 

1 1 

Spacinj 

y 

III 

1/8 

1 1 1 

1/4 

1 1 1 

1/2 

1 1 

1 

Figure  9:  Example  Floating  Point  Model 


The  real  number  line  and  the  floating  point  numbers  which  approximate  it  for 
b = 2,  t = 3,  emin  = “1?  s-nd  Cmax  = 2 are  depicted  in  Figure  9.  The  set  of  positive 
floating  point  numbers  can  be  easily  counted.  For  each  of  the  Cmax  + 1 — ^min  values 
of  the  exponent  e,  there  are  (6  — 1)6*”^  positive  floating  point  numbers.  For  the 
example  of  Figure  9,  there  are  33  floating  point  numbers;  16  are  positive,  16  are 
negative  and  one  is  zero.  These  numbers  are  not  evenly  spaced.  Numbers  with 
small  exponents  are  near  the  origin  and  closely  spaced,  while  numbers  with  large 
exponents  are  far  from  the  origin  and  widely  spaced.  In  fact,  the  absolute  spacing 
between  consecutive  positive  floating  point  numbers  with  exponent  e is  b^b^~*  which 
depends  on  e.  On  the  other  hand,  the  largest  relative  spacing  for  positive  floating 
point  numbers  with  exponent  e is  6^”^  which  is  independent  of  e.  Here,  the  spacing 
has  been  divided  by  6®,  the  smallest  positive  floating  point  number  with  exponent 
e. 

One  particular  parameter  of  interest,  machine  epsilon,  can  be  derived  from  the 
base  b and  the  precision  t.  The  smallest  floating  point  number  that  can  be  added  to 
1.0  so  that  the  result  is  a floating  point  number  larger  than  1.0  is  called  machine 
epsilon.  If  default  rounding  is  used,  this  number  is  Cm  = since  is  the 
spacing  between  consecutive  positive  floating  point  numbers  with  exponent  0.  For 
the  floating  point  system  illustrated  in  Figure  9,  machine  epsilon  is  1/8. 

The  error  in  representing  a nonzero  real  number,  x,  on  the  computer  can  be 
expressed  in  terms  of  machine  epsilon.  This  error  is  sometimes  called  rounding 
error.  Suppose  that  x falls  in  the  range  of  floating  point  numbers;  that  is,  b^  < 
|x|  < (1  — b~^)b^'^^  where  Cnun  ^ e < Cmax-  Using  default  rounding,  the  maximum 
distance  (or  error)  between  x and  the  closest  floating  point  number  is  6^€in.  The 
relative  distance  between  x and  the  closest  floating  point  number  is  less  than  or  equal 
to  Cm-  If  |x|  < X is  said  to  “hard”  underflow.  The  IEEE  Standard  also 

defines  a “soft”  underflow  by  allowing  the  floating  point  numbers  with  exponent  Cmin 
to  be  denormalized;  that  is,  bi  is  allowed  to  be  zero.  In  this  case  “soft”  underflow 
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occurs  if  |a:|  < If  |a;|  > (1  — x is  said  to  overflow.  Compilers 

often  provide  options  for  dealing  with  underflows  and  overflows.  The  default  is 
typically  to  set  underflows  to  zero  and  to  stop  execution  when  overflows  occur. 

A standard  approach  for  obtaining  the  correct  floating  point  characteristics  of 
a computer  is  to  use  the  three  Fortran  function  subprograms,  IlMACH,  for  integer 
parameters,  RIMACH,  for  real  parameters,  and  DIMACH,  for  double  precision 
parameters.  These  routines  were  originally  written  at  BeU  Labs  by  P.  A.  Fox,  A. 
D.  Hall,  and  N.  L.  Schryer  and  are  currently  widely  available  [25].  Many  scientific 
computer  applications  written  in  Fortran  now  have  these  three  function  subprograms 
as  part  of  their  distribution.  When  such  an  application  is  installed  on  a particular 
computer,  the  installer  must  uncomment  the  data  statements  in  these  routines  that 
correspond  to  the  appropriate  computer  before  compilation. 

Numerical  errors  produced  by  an  algorithm  arise  from  errors  in  representing 
the  data,  errors  in  floating  point  arithmetic,  and  amplification  of  errors  due  to  the 
sensitivity  of  various  parts  of  the  algorithm  to  perturbations.  The  notation  fl()  is 
often  used  to  indicate  the  floating  point  result  of  a machine  computation.  If  2 is  a 
real  number  in  the  range  of  the  floating  point  numbers,  then 

11(2)  = Z{1  + €),€<  €ra 

This  means  that  the  relative  error  in  representing  2 on  the  machine  is  bounded  by 
machine  epsilon. 

Let  X and  y denote  nonzero  floating  point  numbers.  The  IEEE  Standard  requires 
that  fl(ar  + y)  equal  the  actual  result,  x + y,  rounded  to  the  nearest  floating  point 
number,  and  similarly  for  subtraction,  multiplication,  and  division.  In  practice,  this 
is  accomplished  by  using  extended  length  registers  for  the  arithmetic.  The  extra 
bits  in  these  registers  are  called  “guard”  bits.  The  number  of  guard  bits  used  varies 
with  vendor.  It  follows  that 


fl(z±y)  = (a:  ± y)(l  + e),  c < Cm 

(30) 

fl(xy)  = zy(l  + c),  6 < Cm 

(31) 

fl(z/y)  = (x/y)(l  + c),  c < Cm. 

(32) 

Suppose  on  the  other  hand  that  x and  y are  the  result  of  a machine  computation 
whose  exact  answer  is  given  by  the  nonzero  numbers  x and  y.  Let  and  €y  denote 
the  relative  errors  in  approximating  x and  y by  z and  y;  that  is,  x = x(l  + €x)  and 
y = y(l  + €y).  Then  substituting  these  relations  for  x and  y into  the  right  hand  side 
of  equations  (30),  (31),  and  (32)  it  follows  that 


|fl(z  ±y)-{x±  y)| 
|z±y| 

|fl(zy)  - zy| 

\xy\ 

|fl(z/y)  - z/y| 

\^/y\ 


< 

< 

< 


iJ^kxKl  + Cm)  + + C„)  + Cm 

kx|  + |fj/|  + + 1(1  "I*  €x)(l  + €y)\€ja 


1 + ft 


+ 1 


1 + fr 


I I ^ I 

■^'l  + cj'" 


(33) 

(34) 

(35) 
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Equations  (34)  and  (35)  show  that  the  error  in  x and  y is  not  amplified  by  multi- 
plication or  division.  Also,  equation  (33)  shows  that  if  cancellation  does  not  occur 
in  X ± y,  then  again  the  error  in  x and  y is  not  amplified.  But  if  cancellation  does 
occur  in  X ± y and  this  number  is  close  to  zero,  then  the  errors  in  x and  y can  be 
greatly  amplified.  This  condition  is  usually  referred  to  as  “catastrophic  cancella- 
tion.” Whenever  possible,  such  cancellations  should  be  avoided  in  the  construction 
of  an  algorithm  although  it  is  often  difficult  to  anticipate  where  they  will  occur. 

Equations  (33)  and  (34)  can  be  used  to  show  that  the  equivalent  expressions 
a(b  — c)  and  ab  — ac  are  also  equivalent  numerically  (both  expressions  have  the 
same  error  bounds).  In  fire  modeling  this  may  arise  when  computing  enthalpy 
differences.  There  is  no  numerical  advantage  to  computing  an  enthalpy  difference 
using  CpT{mi  — m2)  instead  of  CpTm\  — CpT'th2. 

Figure  10  illustrates  the  steps  involved  in  multiplying  two  numbers  on  the 
computer  using  the  example  floating  point  model  given  in  Figure  9.  For  this 
example, Cx  = -05,  Cy  = .1,  = 1/8  and 

|fl(xy)  - xy|/lxy|  « .1453 
which  is  indeed  smaller  than 

kx|  + kyl  + kxCyl  + 1(1  + fx)(l  + €j,)€m|  ^ .2994  . 

Some  vendors  provide  additional  extended  length  registers  for  storage  of  interme- 
diate results.  In  this  case,  intermediate  results  are  not  rounded  to  working  precision 
before  they  are  used. 

The  issue  of  loss  of  accurate  digits  due  to  cancellation  arises  in  the  choice  of 
unknowns  in  a zone  fire  model.  The  presence  of  a fire  in  a building  causes  a pres- 
sure change  in  each  room.  This  pressure  change  is  generally  small  (a  few  Pascals) 
compared  with  the  ambient  pressure  Pamb  (usually  one  atmosphere  or  about  10^ 
Pascals).  The  pressure  differential  across  a door  or  other  opening  drives  the  transfer 
of  mass  and  energy  between  connected  rooms.  Accurate  mass  and  energy  transfer 
computations  require  accurate  pressure  differentials  as  input.  Here  is  a case  where 
reasonable  accuracy  is  desirable  in  a computation  which  is  subject  to  catastrophic 
cancellation.  Let  a and  b denote  pressure  changes  in  adjoining  rooms.  If  the  total 
pressure  in  each  room  is  chosen  as  an  unknown,  then  the  pressure  differential  will 
be  computed  as  (Pamb  + a)  — (Pamb  + &)•  On  the  other  hand,  if  the  pressure  rise 
in  each  room  is  chosen  as  an  unknown,  the  pressure  differential  wiU  be  computed 
as  a — 6.  Noting  that  generally  |a|,  |6|  -C  Pamb?  examination  of  equation  (33)  shows 
that  this  second  choice  is  numerically  superior  to  the  first. 
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x=1.3  y=3.6  **y=4.68 


01/2-  2 4 


Figure  10:  Sample  Floating  Point  Multiplication 
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